###########################################################################
# Lipps & Voeten
# Voting with Putin:
# Gender, LGBT Rights, and Tacit Support for Russia among Europe's Parliamentarians
# Comparative Political Studies

# Appendix replication file --------------------------------------------------
###########################################################################

rm(list = ls())

packages <- c("fixest", "tidyverse", "corrplot", "psych", "broom", "ggplot2", "ggpubr")

for (package in packages) {
  if (!requireNamespace(package, quietly = TRUE)) {
    install.packages(package, dependencies = TRUE)
  }
}

#Packages
library(fixest)
library(tidyverse)
library(corrplot)
library(psych)
library(broom)
library(ggplot2)
library(ggpubr)

df <- readRDS("PACE_Russia_votes.rds")
all <- readRDS("PACE_allvotes.rds")

#Labels of main variables
dict = c(v2pagender = "Women party leadership", government = "Opposition", Equality="Cultural Liberalism",
         female="Female", v2pawomlab = "Women in labor force", v2palgbt = "LGBT rights",
         Y = "Pro-Russia vote", v2paculsup = "Cultural superiority", v2xpa_popul ="Populism (V-Party)",
         v2pariglef = "Econ left-right (V-Party)", v2xpa_illiberal = "Party illiberalism (V-Party)",
         voteid = "Vote ID", natparty_vpartyname = "National Party", signator="Signed LGBT declaration",
         dove="Dovishness", country="Country", party_clean="Party group", year="Year", rile="Left-Right (pm)", 
         eu_position = "European integration (CHES)", lrgen = "Left-Right (CHES)",   lrecon = "Econ Left-Right (CHES)",
         galtan = "Galtan (CHES)", civlib_laworder = "civil liberties (CHES)",   sociallifestyle = "social lifestyle",
         us = "US power (CHES)", russian_interference = "Russian interference (CHES)", postpussy = "Post Pussy Riot Performance",
         signat.all = "Signed LGBT Declaration", signat.notrussia = "Signed LGBT Declaration, not Russia")

#Post conservative turn
cutoff_date <- as.Date("2012-02-21")
df$postpussy <- ifelse(df$docdate >= cutoff_date, 1, 0)

#Break Equality variable into quantiles
df$Qequal <- cut(df$Equality, breaks = quantile(df$Equality,
                                                probs = seq(0, 1, 0.333),
                                                na.rm = T),
                 include.lowest = T,
                 right = F,
                 labels =  c("Conservative", "Medium", "Liberal"))
table(df$Qequal)


#Replication code----------------------------------------------------------

#Table A1: MPs by party group and self-identified gender
sharefem <- df %>% 
  distinct(MP_ID, .keep_all = T) %>%
  group_by(party_clean) %>%
  summarize(male=n()-sum(female), female=sum(female), share = sum(female)/n(), total=n())

#Figure A1: Correlation between main independent variables
dfcor <- df %>% select(v2pagender, Equality, v2xpa_illiberal, v2xpa_popul, v2pariglef, dove, eu_position, lrecon, galtan)

M = cor(dfcor, use="complete.obs")
rownames(M) <- dict[rownames(M)]
colnames(M) <- dict[colnames(M)]

pdf(file = "appendix_figA1.pdf",width = 7,height = 7,onefile=F)
corrplot(M, method = 'number')
dev.off()

#Figure A2: Vote alignment with the Russian delegation on issues of gender equality and LGBT rights prior to 2014
genissues <- all %>% filter(country!="RU" & year<2014 & tit.gender==1 & RusPaceAgree==0)

#MPs based on their gender supporting the PACE majority or Russian minority position
genissues %>% group_by(country) %>%
  summarise(russup=mean(Y),
            russup.up = russup + 1.96* sd(Y)/sqrt(length(Y)),
            russup.low = russup - 1.96 * sd(Y)/sqrt(length(Y))) -> russup.bycountry

pdf(file = "appendix_figA2.pdf",width = 10,height = 9,onefile=F)
ggplot(russup.bycountry, aes(x=russup, y=reorder(country, russup))) +
  geom_bar(stat="identity", color="black", position=position_dodge()) +
  labs(y="", x="Percentage of MPs who align with Russia in vote") +
  theme_light() +
  theme(axis.title = element_text(size=16),
        axis.text = element_text(size=16))
dev.off()

#Table B1: Regression with voteID^country fixed effects
model7 <- as.formula("Y ~   Equality + v2pagender | voteid^country")
model8 <- as.formula("Y ~   Equality + v2pagender + Populist_Illiberal + v2pariglef + dove| voteid^country")
model9 <- as.formula("Y ~   i(female, Equality, ref=0) + female + Equality + v2pagender + Populist_Illiberal + v2pariglef + dove | voteid^country")
model10 <- as.formula("Y ~   i(signat.all, Equality, ref=0) + signat.all + Equality + v2pagender + Populist_Illiberal + v2pariglef + dove | voteid^country")

M7 = feols(model7, df)
M8 = feols(model8, df)
M9 = feols(model9, df)
M10 = feols(model10, df)

etable(list(M7, M8, M9, M10), digits = "r3", dict=dict, 
       file = "TableAll_robust_revision.tex", replace = TRUE, 
       title = "Regression on all votes, all countries, using more restrictive voteID*country fixed effects")

#Table B2: using only signatories to declarations not targetting Russia
model2a <- as.formula("Y ~   signat.notrussia + Equality + v2pagender  | voteid + natparty_vpartyname")
model3a <- as.formula("Y ~   i(female, signat.notrussia, ref=0) + female + signat.all + Equality + v2pagender | voteid + natparty_vpartyname")
model6a <- as.formula("Y ~   i(signat.notrussia, Equality, ref=0) + signat.notrussia + Equality + v2pagender | voteid + country^year")

M2a = feols(model2a, df)
M3a = feols(model3a, df)
M6a = feols(model6a, df)

etable(list(M2a, M3a, M6a), digits = "r3", dict=dict, 
       file = "TableSignatnotrussia_revision.tex", replace = TRUE, 
       title = "Regression on all votes, all countries, using only pro-LGBT declaration that are not targetting Russia")

#Table B3: votes on Russian membership rights
# Testing the effect of individual gender and ideology in a party fixed effects model
model1 <- as.formula("Y ~   female + Equality + v2pagender | voteid + natparty_vpartyname")
model2 <- as.formula("Y ~   signat.all + Equality + v2pagender  | voteid + natparty_vpartyname")
model3 <- as.formula("Y ~   i(female, signat.all, ref=0) + female + signat.all + Equality + v2pagender | voteid + natparty_vpartyname")

#Testing whether the effect of party level variables is robust across country-year
#and estimating the interaction between party ideology and individual gender/ideology
model4 <- as.formula("Y ~   Equality + v2pagender | voteid + country^year")
model5 <- as.formula("Y ~   i(female, Equality, ref=0) + female + Equality + v2pagender | voteid + country^year")
model6 <- as.formula("Y ~   i(signat.all, Equality, ref=0) + signat.all + Equality + v2pagender | voteid + country^year")


df_mem <- df %>% filter(russiamem==1)

M1m = feols(model1, df_mem)
M2m = feols(model2, df_mem)
M3m = feols(model3, df_mem)
M4m = feols(model4, df_mem)
M5m = feols(model5, df_mem)
M6m = feols(model6, df_mem)

etable(list(M1m, M2m, M3m, M4m, M5m, M6m), digits = "r3", dict=dict, 
       file = "Tablemem_revision.tex", replace = TRUE, 
       title = "Regression on votes concerning Russia’s membership rights in PACE, all countries")

#Table B4: votes on Russian domestic transgressions
df_domest <- df %>% filter(international_handc==0)

M1d = feols(model1, df_domest)
M2d = feols(model2, df_domest)
M3d = feols(model3, df_domest)
M4d = feols(model4, df_domest)
M5d = feols(model5, df_domest)
M6d = feols(model6, df_domest)

etable(list(M1d, M2d, M3d, M4d, M5d, M6d), digits = "r3", dict=dict, 
       file = "TableDomestic_revision.tex", replace = TRUE, 
       title = "Regression on Votes concerning domestic transgressions, all countries")

#Table B5: votes on Russian international aggressions
df_internat <- df %>% filter(international_handc==1 & russiamem==0)

M1i = feols(model1, df_internat)
M2i = feols(model2, df_internat)
M3i = feols(model3, df_internat)
M4i = feols(model4, df_internat)
M5i = feols(model5, df_internat)
M6i = feols(model6, df_internat)

etable(list(M1i, M2i, M3i, M4i, M5i, M6i), digits = "r3", dict=dict, 
       file = "TableInternational_revision.tex", replace = TRUE, 
       title = "Regression on Votes concerning international aggressions, all countries")

#Table B6: votes on physical integrity rights
df_phys <- df %>% filter(physint_handc==1)

M1p = feols(model1, df_phys)
M2p = feols(model2, df_phys)
M3p = feols(model3, df_phys)
M4p = feols(model4, df_phys)
M5p = feols(model5, df_phys)
M6p = feols(model6, df_phys)

etable(list(M1p, M2p, M3p, M4p, M5p, M6p), digits = "r3", dict=dict, 
       file = "TablePhys_revision.tex", replace = TRUE, 
       title = "Regression on physical integrity rights votes, all countries")

#Table B7: votes on military conflict
df_mil <- df %>% filter(milit_handc==1)

M1mil = feols(model1, df_mil)
M2mil = feols(model2, df_mil)
M3mil = feols(model3, df_mil)
M4mil = feols(model4, df_mil)
M5mil = feols(model5, df_mil)
M6mil = feols(model6, df_mil)

etable(list(M1mil, M2mil, M3mil, M4mil, M5mil, M6mil), digits = "r3", dict=dict, 
       file = "TableMil_revision.tex", replace = TRUE, 
       title = "Regression on votes relating to military conflict, all countries")

#Table B8: using CHES party indicators
# Testing the effect of individual gender and ideology in a party fixed effects model, including female leadership and party ideology on equality as these are likely to change over time
model1c <- as.formula("Y ~   female + v2pagender + galtan + eu_position + lrecon | voteid + natparty_vpartyname")
model2c <- as.formula("Y ~   signat.all + v2pagender + galtan + eu_position + lrecon | voteid + natparty_vpartyname")
model3c <- as.formula("Y ~   i(female, signat.all, ref=0) + female + signat.all + v2pagender + galtan + eu_position + lrecon| voteid + natparty_vpartyname")

#Testing whether the effect of party level variables is robust across country-year and estimating the interaction with individual gender and ideology
model4c <- as.formula("Y ~   v2pagender + galtan + eu_position + lrecon | voteid + country^year")
model5c <- as.formula("Y ~   i(female, galtan, ref=0) + female + v2pagender + galtan + eu_position + lrecon | voteid + country^year")
model6c <- as.formula("Y ~   i(signat.all, galtan, ref=0) + signat.all + v2pagender + galtan + eu_position + lrecon | voteid + country^year")

#All votes, all parties, all countries
M1c = feols(model1c, df)
M2c = feols(model2c, df)
M3c = feols(model3c, df)
M4c = feols(model4c, df)
M5c = feols(model5c, df)
M6c = feols(model6c, df)


etable(list(M1c, M2c, M3c, M4c, M5c, M6c), digits = "r3", dict=dict, 
       file = "TableCHES_revision.tex", replace = TRUE, 
       title = "egression on all votes, all countries, using CHES party scores instead of V-Party variables")

#Table B9: Russia is mentioned in document title
df_rus <- df %>% filter(russia==1)

M1t = feols(model1, df_rus)
M2t = feols(model2, df_rus)
M3t = feols(model3, df_rus)
M4t = feols(model4, df_rus)
M5t = feols(model5, df_rus)
M6t = feols(model6, df_rus)

etable(list(M1t, M2t, M3t, M4t, M5t, M6t), digits = "r3", dict=dict, 
       file = "TableRus_revision.tex", replace = TRUE, 
       title = "Regression on votes with an explicit reference to Russia in the document title, all countries")

#Table B10: regressions on subset of countries with low threat
df_low <- df %>% filter(country %in% c("BG", "GR", "HU", "IE", "IT", "LU", "ES", "SI", "PT", "AT", "CH", "LI", "FR", "HR", "DE", "UK", "NL", "BE", "DK", "CY", "AD", "SM", "MT", "MC", "TR", "ME", "IS"))

M1low = feols(model1, df_low)
M2low = feols(model2, df_low)
M3low = feols(model3, df_low)
M4low = feols(model4, df_low)
M5low = feols(model5, df_low)
M6low = feols(model6, df_low)

etable(list(M1low, M2low, M3low, M4low, M5low, M6low), digits = "r3", dict=dict, 
       file = "TableWest_revision.tex", replace = TRUE, 
       title = "Regression on all Votes, US allies for whom Russia is not a high threat")
subsampresults1 <- tidy(M3low) %>% mutate(model="Russia not high threat")

#Table B11: regressions on subset of countries with high threat
df_high <- df %>% filter(country %in% c("CZ", "NO", "RO", "SK", "SE", "FI", "EE", "LT", "LV", "PL", "UA", "GE", "MD", "AL", "MK"))

M1h = feols(model1, df_high)
M2h = feols(model2, df_high)
M3h = feols(model3, df_high)
M4h = feols(model4, df_high)
M5h = feols(model5, df_high)
M6h = feols(model6, df_high)

etable(list(M1h, M2h, M3h, M4h, M5h, M6h), digits = "r3", dict=dict, 
       file = "TableEastNATO_revision.tex", replace = TRUE, 
       title = "Regression all votes, US allies for whom Russia is high threat")
subsampresults2 <- tidy(M3h) %>% mutate(model="Russia high threat")

#Table B12: regressions on subset of countries allied with Russia
df_ally <- df %>% filter(country %in% c("RU", "RS", "AM", "AZ", "BA"))

M1al = feols(model1, df_ally)
M2al = feols(model2, df_ally)
M3al = feols(model3, df_ally)
M4al = feols(model4, df_ally)
M5al = feols(model5, df_ally)
M6al = feols(model6, df_ally)

etable(list(M1al, M2al, M3al, M4al, M5al, M6al), digits = "r3", dict=dict, 
       file = "TableEast_revision.tex", replace = TRUE, 
       title = "Regression on all Votes, Russian Allies")

#Figure B1: Interaction of female and signatory in country subsamples
subsampresults <- tidy(M3al) %>%
  mutate(model="Russian ally") %>%
  bind_rows(subsampresults1, subsampresults2) %>%
  filter(str_detect(term, "signat.all") | str_detect(term, "female"))

pdf(file="appendix_figB1.pdf", width=9, height=6)
ggplot(subsampresults, aes(x = term, y = estimate, ymin=estimate-(1.96*std.error), ymax=estimate+(1.96*std.error), shape=model, color=model)) +
  geom_errorbar(width=0.2, size=1, position="dodge") +
  geom_point(position = position_dodge2(width=0.2, preserve = "total"), stat="identity", size=4) +
  geom_hline(yintercept=0) +
  scale_x_discrete(labels = c("Female", "Female * Signatory", "Signatory")) +
  scale_color_manual(values=c("black", "grey48", "grey80"), name="Country sample")+
  scale_shape_discrete(name="Country sample") +
  labs(x = "", y = "Estimate and 95% Conf. Int.", title = "") +
  theme_bw() +
  theme(
    axis.text.x = element_text(size = 16),
    axis.text.y = element_text(size = 16),
    legend.position = "bottom",
    legend.text = element_text(size = 16),
    panel.grid.major = element_line(colour = 'grey', linetype = 'dashed'),
    title = element_text(size = 16, vjust = 0.5))
dev.off()

#Figure B2: Interaction between female and country
Mc=feols(Y~female*country |voteid, df)
coef_df <- as.data.frame(coef(Mc))
coef_names <- rownames(coef_df)

# main effect of female
main_female_effect <- coef_df["female", 1]

# Extract the interaction terms (female:country)
interaction_terms <- coef_names[grepl("female:country", coef_names)]
interaction_effects <- coef_df[interaction_terms, , drop = FALSE]

# main effects for countries
country_terms <- coef_names[grepl("^country", coef_names) & !grepl("female:country", coef_names)]
country_effects <- coef_df[country_terms, , drop = FALSE]

# data frame for plotting
plot_data <- data.frame(
  country = gsub("female:country", "", interaction_terms), 
  interaction_effect = interaction_effects[, 1],
  se_interaction = sqrt(diag(vcov(Mc))[interaction_terms]))

# Add main effect of female to interaction effects
plot_data <- plot_data %>%
  mutate(main_female_effect = main_female_effect,
    total_effect = main_female_effect + interaction_effect,
    se_main = sqrt(diag(vcov(Mc))["female"]),
    lower = total_effect - 1.96 * sqrt(se_main^2 + se_interaction^2),
    upper = total_effect + 1.96 * sqrt(se_main^2 + se_interaction^2))

# Sort countries by total effect for better visualization
plot_data <- plot_data %>%
  arrange(total_effect)

pdf(file="appendix_figB2.pdf", width=8, height=9)
ggplot(plot_data, aes(x = reorder(country, total_effect), y = total_effect)) +
  geom_point(size=1.5) +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.3, linewidth=1) +
  geom_vline(xintercept=0, linetype="dashed") +
  coord_flip() +
  theme_minimal() +
  labs(x = "", y = "Total Effect of Female (with 95% CI)") +
  theme(axis.text = element_text(size = 16),
        axis.title = element_text(size = 16),
        panel.grid.major = element_line(colour = 'grey80', linetype = 'dashed'))
dev.off()

#Figure B3: Female party leadership pre- and post-Pussy riot
modelpussy2 <- as.formula("Y ~  v2pagender*postpussy | voteid + natparty_vpartyname")
M1pussy2 = feols(modelpussy2, df)
M2pussy2 = feols(modelpussy2, df_domest)
M3pussy2 = feols(modelpussy2, df_internat)

pdf(file = "appendix_figB3.pdf", width=8, height=7, onefile = F)
coefplot(list(M1pussy2, M2pussy2, M3pussy2),
         ci.col=c("black", "grey48", "grey80"), pt.col= c("black", "grey48", "grey80"), ci.width=0.04, cex= 2, lwd=2.75,
         dict=c("v2pagender"="Women party leadership", "v2pagender:postpussy"="Women party leadership x Post-Pussy Riot"), main="")
legend("topright", col = c("black", "grey48", "grey80"), pch = c(16,17,15), lwd = 1.25, legend = c("All Votes", "International Votes","Domestic Votes"), title = "Type of Votes")
dev.off()

#Table B13: the effect of female on signing a LGBT rights declaration
M1signat = feols(signat.all ~ female | year^natparty_vpartyname , df) 
M2signat = feols(signat.all ~ i(female, Qequal) | year^country , df)
M3signat = feols(signat.notrussia ~ female | year^natparty_vpartyname, df)
M4signat = feols(signat.notrussia ~ i(female, Qequal) | year^country , df)

etable(list(M1signat, M3signat, M2signat, M4signat), digits = "r3", dict=dict, 
       file = "TableDVsignat_revision.tex", replace = TRUE, 
       title = "The effect of gender on signing a declaration condemning the violation of LGBT rights")

#Figure B4: the effect of female on signing a LGBT rights declaration
signatresults1 <- tidy(M2signat) %>% 
  mutate(DV = "All Declarations")
signatresults2 <- tidy(M4signat) %>%
  mutate(DV= "Declarations not on Russia") %>%
  bind_rows(signatresults1) %>%
  mutate(liberalism_level = case_when(
    str_detect(term, "Qequal::Conservative") ~ "Conservative",
    str_detect(term, "Qequal::Medium") ~ "Moderate",
    str_detect(term, "Qequal::Liberal") ~ "Liberal")) %>%
  mutate(liberalism_level = factor(liberalism_level, levels = c("Conservative", "Moderate", "Liberal")),
         female = case_when(str_detect(term, "female::0") ~ 0,
                            TRUE ~ 1))

pdf(file = "appendix_figB4.pdf",width = 10,height = 7,onefile=F)
ggplot(signatresults2, aes(x = liberalism_level, y = estimate, color = as.factor(female), ymin=estimate-(1.96*std.error), ymax=estimate+(1.96*std.error))) +
  geom_errorbar(width=0.3, size=1, position="dodge") +
  geom_point(position = position_dodge2(width=0.3, preserve = "total"), stat="identity", size=2) +
  geom_hline(yintercept=0, color="black") +
  scale_color_manual(labels= c("Male", "Female"), values=c("deeppink4", "deepskyblue4"), name="") +
  labs(x = "Party Cultural Liberalism", y = "Coefficient Estimate", title="Effect on signing an pro-LGBT declaration") +
  facet_wrap(~DV)+
  theme_bw() +
  theme(axis.text.x = element_text(size=16), axis.text.y=element_text(size=16), 
        legend.position = "bottom", 
        legend.text = element_text(size=16),
        panel.grid = element_line(colour = 'grey', linetype = 'dashed'), 
        title = element_text(size=16, vjust=0.5),
        strip.text = element_text(size=16))
dev.off()

#the results with coarsened exact matching can be found in the separate script "replication_matching.R"